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In this paper we discuss the anatomy of frequency-domain gravitational-wave signals from non-precessing 
black-hole coalescences with the goal of constructing accurate phenomenological waveform models. We first 
present new numerical-relativity simulations for mass ratios up to 18, including spins. From a comparison 
of different post-Newtonian approximants with numerical-relativity data we select the uncalibrated SEOBNRv2 
model as the most appropriate for the purpose of constructing hybrid post-Newtonian/numerical-relativity wave¬ 
forms, and we discuss how we prepare time-domain and frequency-domain hybrid data sets. We then use our 
data together with results in the literature to calibrate simple explicit expressions for the final spin and radi¬ 
ated energy. Equipped with our prediction for the final state we then develop a simple and accurate merger- 
ringdown-model based on modified Lorentzians in the gravitational wave amplitude and phase, and we discuss 
a simple method to represent the low frequency signal augmenting the TaylorE2 post-Newtonian approximant 
with terms corresponding to higher orders in the post-Newtonian expansion. We finally discuss different options 
for modelling the small intermediate frequency regime between inspiral and merger-ringdown. A complete 
phenomenological model based on the present work is presented in a companion paper Q 
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I. INTRODUCTION 

Our goal is to develop an accurate and simple description 
of the gravitational-wave (GW) signal of a compact binary co¬ 
alescence in the form of explicit expressions in the frequency 
domain, which are particularly convenient for GW data anal¬ 
ysis. Optimal methods for the detection and accurate identi¬ 
fication of events are based on filtering the data stream with 
a template bank of physically correct waveforms. Extract¬ 
ing the maximal amount of information from the data thus 
relies critically on the availability of accurate waveform mod¬ 
els, such as the ones we discuss here. For small masses, such 
as for neutron-star binaries, the post-Newtonian (PN) descrip¬ 
tion El is considered sufficient for searches for GW events in 
the advanced detector era usm. However, as the binary mass 
increases and the actual merger comes into band, the PN ap¬ 
proximation breaks down and non-perturbative solutions for 
the Einstein equations are required. Eor the simple case of 
non-spinning black holes (BHs), a comparison of different PN 
results suggests that such non-perturbative information is re¬ 
quired for the efficient detection of binaries with a total mass 
larger than 12 solar masses i). In this work we will treat the 
case of coalescing BHs, or neutron stars spiralling into suffi¬ 
ciently heavy BHs such that equation of state effects can be 
neglected. The final state of the coalescence, a single per¬ 
turbed Kerr BH, can also be treated with linear perturbation 
theory. This allows us to compute the complex frequencies of 
damped sinusoid spheroidal harmonic “quasinormal” modes, 
but not their amplitudes or relative phase. 

In the present paper we focus on understanding the 
anatomy of frequency-domain waveforms and laying the 
groundwork for accurate frequency-domain inspiral-merger- 


ringdown (IMR) waveform models, including a discussion of 
how we combine information from numerical relativity and 
perturbative approaches. In a second paper [I] we describe 
in detail the construction of a specific model, “PhenomD”, 
which we have implemented in the LIGO Algorithms Library 
(LAE) 0 , the key software infrastructure for GW data analy¬ 
sis from interferometric ground based detectors, and we com¬ 
pare with other waveform models implemented in LAE, in 
particular SEOBNRv2 0. While the main motivation is to 
provide a practical tool for GW astronomy that we hope to 
be useful for searches or parameter estimation, we stress that 
gaining a better understanding of waveforms is by itself valu¬ 
able for GW physics. 

Astrophysical models of compact binary merger event rates 
lead to the expectation of a first direct detection of GWs within 
the next five years cm, as a new generation of interferomet¬ 
ric detectors is approaching design sensitivity EHH, expand¬ 
ing the volume of the universe observable in the GW window 
by roughly a factor of 1000 over previous detectors. This in¬ 
crease is expected to be sufficient to exhaust the current large 
uncertainties in event rates, which are connected to the large 
uncertainties about the populations of compact binaries in the 
universe. Our lack of understanding of binary populations 
calls for being prepared to cover as much as possible of the 
physical binary black hole (BBH) parameter space, and the 
waveform model we discuss has been calibrated to the largest 
region of the non-precessing parameter space to date, covering 
mass ratios up to 18. 

The first successful calculation of the GW signal during 
roughly the last orbit, merger and subsequent quasi-notmal- 
mode (QNM) ringdown from a numerical solution of the Ein¬ 
stein equations, without perturbative approximations, dates 
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back only one decade im. Since then, synthesizing models of 
the complete coalescence waveform from PN results, numer¬ 
ical relativity (NR) and BH perturbation theory has become 
a key goal of GW source modeling. One approach has been 
to model gravitational waveforms in the time domain, by cali¬ 
brating an effective-one-body (EOB) resummation Cilia of 
PN results to NR simulations ll^ [T6l420l . This approach pro¬ 
vides very accurate waveforms, but still requires the solution 
of a system of ordinary differential equations for each case and 
is computationally expensive, in particular for applications in 
parameter estimation and model selection within a Bayesian 
framework. A route to speed up the evaluation of EOB models 
by means of reduced order techniques is discussed in GIIISI. 

Our previous work has followed an alternative approach, 
where we have modeled the wave signal with explicit expres¬ 
sions in the frequency domain. This approach offers simplic¬ 
ity and evaluation speed, however previous models showed 
relatively poor extrapolation beyond their calibration region 
and limited accuracy for parameter estimation. In this work 
we redesign our phenomenological approach based on a more 
detailed study of the anatomy of waveforms in the frequency 
domain and calibrate the model with a larger and more ac¬ 
curate set of numerical waveforms up to a mass ratio of 
18. As for previous non-precessing phenomenological wave¬ 
form models Il24l427l . this work is restricted to the dominant 
{ - \m\ - 2 spherical harmonic mode, leaving the general 
aligned spin case, higher harmonics and precession for future 
work. Erom the point of view of GW data analysis, the model 
we construct here can directly replace our previous PhenomB 
and PhenomC waveform models with similar computational 
cost but an accuracy that for a large part of the parameter space 
further improves upon current EOB-NR models. A route to¬ 
ward extending such models to precession has been described 
in 1281 . 


In the next section, we will discuss our “input” waveforms; 
NR waveforms, inspiral waveforms computed in the EOB de¬ 
scription, and other PN approximants, and how we construct 
complete time domain hybrid waveforms. Our approach does 
in fact allow us to model the waveforms in pieces, e.g. to 
calibrate the inspiral and merger-ringdown models separately, 
without ever constructing a hybrid first. However, the con¬ 
struction of complete time domain hybrids is convenient for 
comparisons, and to develop our modeling approach without 
artificial restrictions due to the length of available waveforms. 
At the hybridisation stage we will compare the agreement 
of different PN approximants with our numerical data, and 
choose to hybridize with “uncalibrated” SEOBNRv2 wave¬ 
forms, i.e., where we have removed all NR calibrations. 


In Sec. Ill we will refocus from the inspiral to the mod¬ 
elling of the final state and present fits for the final spin and 
radiated energy as functions of mass ratio and spin. These fits 
will provide crucial information regarding the modelling of 
the ringdown. With the final state information and a complete 
time domain hybrid in hand, we will turn to the study of the 
anatomy of waveforms in the Eourier domain which will moti¬ 
vate our plan for construction of a waveform model, presented 
in Sec. IV We conclude with a summary in Sec. |V] 


II. INPUT WAVEFORMS AND HYBRID CONSTRUCTION 
A. Waveform notation and conventions 


Our work concerns the ( - \m\ - 2 spherical harmonic 
modes of the complex GW strain h 22 , which we decompose 
into an amplitude A and phase (j), which depend on time f, and 
the intrinsic parameters of the source E, 

( 2 . 1 ) 


The intrinsic parameters S consist of the dimensionless pro¬ 
jections of the BH spins § 1^2 in the direction of the orbital 
angular momentum L, and the masses mi_ 2 , where 


A; = 


^i-L 


( 2 . 2 ) 


and we define the mass ratio q - oti/ot 2 ^ 1. total mass M - 
nil +'«2 and symmetric mass ratio 77 = 17111112 !M^. Our binaries 
are non-precessing, and as usual we choose our coordinates 
and spherical harmonic modes consistent with the equatorial 
symmetry exhibited by such spacetimes, such that 


h22 = /J2,_2- (2-3) 

We also assume negligible eccentricity, in which case the fre¬ 
quency (xi{t) - dip I dt is a monotonic function of f. Resid¬ 
ual eccentricity in numerical waveforms and potentially other 
numerical artefacts will cause oscillations, however, pro¬ 
vided these are small enough, monotonicity is preserved until 
merger, and we can define the inverse function t(a)). 

We define the Eourier transform Ii{(jj) of a complex function 
h(t) as 


h(f)^ I dt. 


(2.4) 


consistent with the conventions adopted in the LIGO Algo¬ 
rithms Library Q. We will also use the radial frequency 
oj - 2nf, in particular in quantitative statements, to facili¬ 
tate comparisons with data analysis considerations. With our 
convention of the Eourier transform, time derivatives are con¬ 
verted to multiplication in the Eourier domain by i2jif. We 
recall, that as a consequence of time derivatives being con¬ 
verted to multiplication in Fourier space, the standard conver¬ 
sion between GW strain and the Newman-Penrose scalar 1 / 14 , 
where 


d^h{t) 

dfi 




(2.5) 


only affects the Fourier domain amplitude, but not the phase, 
up to a jump of tt, and apart from possible effects specific to 
the numerical algorithm used to carry out this conversion. 

For NR data or hybrid waveforms we need to Fourier trans¬ 
form discrete time series hr, then the above convention implies 


n-l 

r=0 


( 2 . 6 ) 
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where Af is the uniform time step of the time series. 

The low frequency limit for the Fourier domain waveform 
can be understood from PN theory. Analytically Fourier- 
transforming the TaylorT2 approximant with the stationary- 
phase-approximation yields the TaylorF2 approximant, an ex¬ 
plicit expression of the Fourier domain strain in the Frequency 
domain. To leading order the amplitude diverges at low fre¬ 
quency as and the phase as see the discussion in 

Secs. IIVAII and IIVA2l below. 

For the modelling of IMR waveforms we also need to un¬ 
derstand the falloff properties of the waveform at high fre¬ 
quencies. For smooth square integrable functions of time it is 
well known (see e.g. mi) that 

kf)^0{\fV^) as l/l^cx, for all M, (2.7) 

thus the Fourier transform then falls off faster than any power 
law at high frequency. If a function h{t) only has p continuous 
derivatives in for some p > 0 and a pth derivative in of 
bounded variation, then mi 

kf) ^ 0(\fr-^) as l/l^cx,. (2.8) 

While the Fourier domain strain diverges as at low fre¬ 
quencies, and is thus not square integrable for idealized in¬ 
finitely long signals (in contrast to finite astrophysical sig¬ 
nals), the time derivative of the strain (the news function) is 
square integrable (being proportional to / at low frequen¬ 
cies), and the faster-than-power-law falloff for the strain can 
be inferred from the scaling of the Fourier amplitude when 
taking time derivatives. 

B. Numerical relativity waveforms 

Our numerical waveforms have been taken from two inde¬ 
pendent data sets, the publicly available SXS catalogue EqI 
(computed by the SpEC code BTHJTII I. and from a set of 
waveforms that have recently been constructed with the BAM 
code 1381 mi . and which are listed in tables||]|II] Both the SXS 
and BAM data sets are divided into waveforms we use to cal¬ 
ibrate the PhenomD model, and a larger set we use for model 
verification. The key data sets that determine the calibration 
range of our waveform model are the BAM waveforms for a 
range of spins at mass ratio 1:18, high-spin BAM data at mass 
ratios 4 and 8, and equal mass SXS data sets at very high spins 
of -0.95 and -hO.98. 

SpEC is a multi-domain pseudo-spectral code that uses ex¬ 
cision techniques to remove the BH interiors (and, most im¬ 
portantly, each BH singularity) from the computational do¬ 
main. Eor an in-depth explanation see Refs. l30l l40l 1411 . 
The SpEC code evolves the Einstein field equations using the 
Generalised Harmonic coordinate formulation. This code has 
been used to produce the longest simulations to date (175 or¬ 
bits / 350 GW cycles for a mass-ratio 1:7 configuration BTIO . 
and 201 waveforms have been made publicly available as of 
June 2015 imil : see Ref. ll3^ for details. 

A key set of waveforms in the SXS catalogue are those for 
highly spinning equal-mass BH binaries. Most codes in the 
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T2 

D/M 

p,IM 

-PrIM 

(xlO-'*) 

Mali 

Ngw 

2 

0.50 

0.50 

12.2 

0.07277 

3.468 

0.04128 

26.8 

2 

0.75 

0.75 

11.1 

0.07591 

3.404 

0.04689 

23.5 

3 

-0.50 

-0.50 

12.8 

0.06309 

2.870 

0.03958 

19.9 

4 

-0.75 

-0.75 

12.1 

0.05694 

1.858 

0.04283 

15.6 

4 

-0.50 

-0.50 

11.9 

0.05669 

2.694 

0.04352 

18.1 

4 

-0.25 

-0.25 

11.5 

0.05694 

2.784 

0.04528 

18.9 

4 

0.00 

0.00 

11.2 

0.05710 

2.829 

0.04673 

21.1 

4 

0.25 

0.25 

11.0 

0.05681 

2.962 

0.04785 

23.0 

4 

0.50 

0.50 

10.8 

0.05648 

4.493 

0.04870 

26.1 

4 

0.75 

0.75 

10.7 

0.05604 

1.136 

0.04931 

30.0 

8 

-0.85 

-0.85 

10.0 

0.04135 

2.6142 

0.06486 

8.2 

8 

0.80 

0.00 

8.00 

0.04146 

3.009 

0.07266 

23.3 

8 

0.85 

0.85 

6.50 

0.04703 

8.0706 

0.10951 

15.7 

10 

0.00 

0.00 

8.39 

0.03670 

1.685 

0.06942 

13.4 

18 

-0.80 

0.00 

10.0 

0.02080 

0.6589 

0.05614 

14.3 

18 

-0.40 

0.00 

9.00 

0.02176 

0.7842 

0.06395 

15.1 

18 

0.00 

0.00 

7.58 

0.02379 

1.261 

0.07932 

13.0 

18 

0.40 

0.00 

7.43 

0.02300 

1.161 

0.08052 

23.2 


TABLE I: Details of new BAM simulations. For each mass ratio q 
and choice of spins andxi, the initial binary separation is D/M, 
and the tangential and orbital momenta are (p,,p,.)IM. The initial 
GW frequency is Mw, , and the simulation completes Agw orbits be¬ 
fore merger. 


g 

Ti 

Xi eccentricity Mf 
(xlO-3) 

a; 

2 

0.50 

0.50 

1.2 

0.945 

0.805 

2 

0.75 

0.75 

4.4 

0.930 

0.889 

3 

-0.50 

-0.50 

1.0 

0.9779 

0.300 

4 

-0.75 

-0.75 

0.8 

0.9846 

0.049 

4 

-0.50 

-0.50 

1.0 

0.9831 

0.194 

4 

-0.25 

-0.25 

1.0 

0.981 

0.334 

4 

0.00 

0.00 

1.4 

0.978 

0.472 

4 

0.25 

0.25 

2.4 

0.9738 

0.607 

4 

0.5 

0.5 

3.6 

0.9674 

0.738 

4 

0.75 

0.75 

4.0 

0.9573 

0.863 

8 

-0.85 

-0.85 

0.5 

0.9931 

-0.320 

8 

0.80 

0.00 

4.9 

0.977 

0.860 

8 

0.85 

0.85 

9.1 

0.9746 

0.895 

10 

0.00 

0.00 

0.8 

0.9918 

0.255 

18 

-0.80 

0.00 

0.5 

0.9966 

-0.531 

18 

-0.40 

0.00 

0.5 

0.9966 

-0.188 

18 

0.00 

0.00 

1.3 

0.9959 

0.163 

18 

0.40 

0.00 

1.8 

0.9943 

0.505 


TABLE II: Eccentricity measured from the orbital frequency, final 
Kerr parameter and final mass of the new BAM simulations. 


field (including BAM and SpEC) primarily use conformally 
flat initial data, which limits the BH spins to a jm ~ 0.92 m- 
l45l . There has been some work in constructing high-spin data 
for puncture codes like BAM ll46l l47l . but to date no wave¬ 
forms for long-term quasi-circular inspirals have been pub¬ 
lished. Eor the SpEC code, however, high-spin data have been 
produced and used in production simulations, based on a su¬ 
perposition of Kerr-Schild data ll^ l44l l48l . This technique 
has made possible a number of equal-mass binary simulations 
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Code/ID ;t'i ;(-2 Mf Of M/rd M/hyb 

SXS:BBH:0001 0.00 0.00 0.9516 0.6865 0.0881 0.00398 

SXS:BBH:0156 -0.95 -0.95 0.9681 0.3757 0.0713 0.00522 

SXS:BBH:0172 0.98 0.98 0.8892 0.9470 0.1328 0.00497 


TABLE III: High spin equal mass waveforms from the SXS cata¬ 
logues used in this paper to illustrate waveform anatomy. For each 
configuration we list the SXS catalogue number, the spins and^ 2 - 
The final BH has mass Mf and dimensionless spin a/, and the QNM 
ringdown signal has frequency The frequency M/^yb marks 

the midpoint of the transition region between SEOBv2 inspiral and 
NR data. 


with high spins, and these have been used here both in the cal¬ 
ibration of our model, and in tests of its accuracy and fidelity. 

The SXS waveforms that we use all cover 12 orbits or more 
of inspiral before merger, and the waveforms are extrapolated 
to null infinity at third order as discussed in ll49l . In this paper 
we will use 3 waveforms as listed in Table [111] from the SXS 
catalogue to illustrate waveform anatomy, and for comparison 
with different PN approximants. 

The BAM code Il38l solves the 3 - 1-1 decomposed Ein¬ 
stein evolution equations using the ;^f-variant of the moving- 
punctures version of the BSSN ifSOlfSTIl formulation. Spatial 
derivatives are sixth-order accurate in the bulk [391. Kreiss- 
Oliger dissipation terms converge at fifth order, and a fourth- 
order Runge-Kutta algorithm is used for the time evolution. 
BBH puncture initial data ll5^ l5Tl are calculated with the 
pseudo-spectral elliptic solver described in Ref. IH). The 
GWs are calculated using the Newman-Penrose scalar T' 4 . For 
further details see Ref ITSl . The GWs are extracted at a finite 
distance, typically ~ 100 M from the source. 

For this work, new simulations were performed at mass ra¬ 
tios 2, 3, 4, 8 and 18. We also repeat the nonspinning 1:10 
simulation reported in Refs. ElEgl, using the same parame¬ 
ters as given in that work. 

For the simulations at mass ratios < 4, the eccentricity 
was reduced to e < 10 “^ using the procedure described in 
Ref El, which is an iterative procedure based on compar¬ 
isons against PN and FOB inspirals. At mass ratios 8 and 18, 
some of the simulations were too short to reliably use this pro¬ 
cedure. However, for large mass ratios the time scale of sig¬ 
nificant early time perturbations of the orbital quantities due to 
gauge transients is smaller than for comparable masses, which 
allows to easily relate observed residual eccentricity in the or¬ 
bital motion to an excess or deficiency in the value of the ini¬ 
tial tangential component of the black hole momentum, which 
dominates the error with respect to the true quasi-circular in¬ 
spiral parameters. Consequently, a straightforward iteration 
“by hand” can be performed rather easily, starting with pa¬ 
rameters calculated from PN inspirals. 

The length of our NR waveforms, summarized in tables |I] 
[H] was chosen partly to sensibly connect them to inspiral ap¬ 
proximants. PN inspiral waveforms become less accurate at 
higher frequencies, and we need to decide at what frequency 
we need to be able to switch to NR information in order to 
produce a sufficiently accurate model. Previous work sug¬ 
gests that ~5-10 orbits will be sufficient to produce models 


that meet the needs for GW detections with advanced GW de¬ 
tectors ll58l I 59 I . if we use standard PN approximants for the 
inspiral. Those works took the largest differences between the 
highest-order PN waveforms available at that time as a con¬ 
servative estimate of their overall error. Here, however, we 
push the boundaries of the NR-accessible parameter space to 
higher mass ratios and spins, and we aim to produce a wave¬ 
form model suitable not only for detection but also parameter- 
estimation purposes in the advanced detector era. 


If we were to follow the stringent accuracy requirements 
for combining (any) PN approximant with NR data, we would 
need to simulate hundreds of NR orbits (60H63l, which is 
prohibitively expensive, especially given the parameter space 
that we are covering. At mass ratio 18, medium-resolution 
simulations of only 8 orbits required 800000 CPU hours and 
high-resolution simulations 1.3 million CPU hours. Previous 
length-requirement studies, however, conservatively treated 
all PN approximants as equally (in)accurate, and it is worth in¬ 
vestigating whether an individual inspiral description is more 
accurate than others, thereby allowing a matching with NR 
data at higher frequencies. While this question is difficult to 
answer exhaustively without extremely long NR waveforms, 
we discuss in Sec. IID that the dramatically better agreement 
between NR and members of the FOB family complements 
the evidence that FOB is a more accurate inspiral description 
Ensnia, hence we can connect it reliably to our reason¬ 
ably short high-mass-ratio simulations. 


We must also ensure that our simulations are sufficiently 
accurate. In the past we have performed detailed convergence 
tests of our code ESI, which provide strong evidence that the 
numerics are correct. Since performing the simulations re¬ 
ported in Ref. (66), we have found that the small number of 
AMR buffer zones (six) limited the accuracy in long unequal- 
mass simulations. This is the cause of the disagreement be¬ 
tween the BAM and SpEC nonspinning q = 4 simulations 
that were discussed in Sec. IV.A of Ref. m- In our current 
simulations we use a minimum of 16 buffer zones, which is 
sufficient to buffer mesh refinement boundary effects during 
a fourth-order Runge-Kutta time step of a fine grid for sixth 
order finite differencing and compatible Kreiss-Oliger dissi¬ 
pation ll67l . 


Even for the lowest-cost configuration (an equal-mass, non 
spinning binary), the computational cost prohibits a standard 
factor-of-two convergence test, and for many configurations 
there have been no clean demonstrations of convergence for 
any BBH code. Our confidence in NR results in general 
is based on the strong agreement between simulations from 
many independent codes, e.g., Refs. Il68l470l . The results of 
a convergence test for a series of simulations for the q' = 18 
nonspinning case are shown in Fig. [T] Due most likely to 
poorly resolved initial gauge dynamics and an initial burst of 
unphysical GW content (“junk radiation”), we have not been 
able to show convergence of the phase aligning at the begin¬ 
ning of the simulation, but we see reasonable sixth order con¬ 
vergence when aligning the i - m - 2 GW signal near the 
merger at a frequency of Mojcw - 0.2. Fig. shows results 
from grid configurations with (96,120,144) points in the in¬ 
nermost mesh refinement boxes, which dominate phasing er- 
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t/M 

FIG. 1: Convergence and Richardson-extrapolated error for non¬ 
spinning BAM simulations at ^ = 18 at resolutions (96,120,144) 
as described in the main text. 


5 0.0120 
S 0.0115 
^ 0.0110 
g 0.0105 
1 0.0100 
I 0.0095 
^ 0.0090 

§ 2 4 6 8 10 12 

convergence order 

FIG. 2: RMS difference between different resolutions as defined in 
Eq. \2.9) for non-spinning BAM simulations at ^ = 18 at resolutions 
(96,120,144) as described in the main text. Rescaling for 6th order 
convergence as shown in Fig.[^is most consistent with our data. 

rors and computational cost. Differences are scaled to sixth 
order convergence and we also show the error obtained as the 
difference of the highest resolution (144) and the Richardson- 
extrapolated result. Secular dephasing is very small, below 
O.QArad, and errors are dominated by oscillations and noise. 
Due to the oscillations and noise it is hard to decide by eye 
which convergence order is most consistent with the data, so 
we considered the RMS deviation 


a spinning q — IS case is practically indistinguishable for dif¬ 
ferent resolutions. We have therefore not carried out a higher 
resolution run for the spinning case to save another 2 mil¬ 
lion CPU hours, and we use the 96 points resolution or better 
for the q' = 18 production runs. 

Our NR waveforms only model the waveform from an as- 
trophysical inspiral after some time Iq, when the influence of 
initial transients can be neglected. When processing a large 
number of waveforms, it is desirable to automatize the tedious 
procedure to find appropriate choices of fo, and we thus need 
to understand initial transients in the waveform well enough 
to exclude them with a general algorithm. Two types of tran¬ 
sients are of particular interest; the unphysical radiation con¬ 
tent in the initial data, often referred to as “junk radiation”, 
and a transient component resulting from our conversion of 
to strain. In order to estimate the time on which junk radiation 
degrades the quality of the wave signal, we first determine the 
maximum of the amplitude of the junk radiation tj for the NR 
1/^4 waveforms, and then add a multiple of the QNM damping 
period Tro - Ij fdamp corresponding to the two initial BHs; 

ATj^tj + ajTRD- ( 2 . 10 ) 

In order to convert 1/^4 waveforms to strain we use the fixed 
frequency integration algorithm lITTII . which depends on the 
choice of cutoff frequency /q. We choose this frequency as 
three quarters of the initial orbital frequency, which for the 
SXS data sets is determined from the included metadata, and 
for the BAM data sets is computed from the PN initial data pa¬ 
rameters for the BH momenta. This value of the initial orbital 
frequency is cross-checked from the actual numerical wave¬ 
form data by extrapolating the frequency evolution of the nu¬ 
merical waveform to the initial time by fitting to a fourth order 
polynomial in a time window starting after a time Tj. As can 
be expected on dimensional grounds, this procedure causes 
transients in the strain that persist for a time scale of the order 
of I// 0 . We thus define 

AT, = crj/o, (2.11) 

and discard numerical data before a time fo, where 

fo = max(ATh,ATj). (2.12) 

For short waveforms of less than 3000M time to peak, we 
choose cTj = l,cr/, = 2.5, and for longer waveforms we 
choose (Tj - 3, cTf, = 1.5. 



RMS^ = f ' ((Xi - X2)/C(n) - (X2 - X2)f dt, 

h -1\ Jf, 

(2.9) 

where the A, are the phases at resolutions (96,120,144) 
aligned at Mtocw - 0.2, C is the rescaling factor correspond¬ 
ing to convergence order n, f j and f 2 are the minimal and max¬ 
imal times shown in Fig. [T] and the unresolved pulse of noise 
between times -500 and -400 is set to zero. As can be seen 
in Fig.J^ the best-fitting convergence order is 6 . In addition, 
in Figj^we show the that the orbital motion for the nonspin¬ 
ning convergence series and a low and medium resolution for 


C. Post-Newtonian and effective-one-body descriptions 

Our work uses PN waveforms in two ways: First, to 
construct “hybrid waveforms” that represent the complete 
inspiral-merger-ringdown of a binary system by appropriately 
gluing together PN and NR waveforms. Without restricting 
the generality of our approach we will carry out this pro¬ 
cedure in the time domain (for an alternative frequency do¬ 
main description see Gzl). Second, we will tune a frequency- 
domain PN approximant to the hybrid waveforms by calibrat¬ 
ing higher order PN terms (which have not yet been computed 
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FIG. 3: Left: Orbital motion of 3 highest resolutions (96, 120, 144) for non-spinning ijr = 18 case, aligned at aligned at MoJorb = 0.2. Right: 
Orbital motion of 2 highest resolutions (96, 120) for ^ = 18, = 0.4,0 case, aligned at MoJorb = 0.1. 


in PN theory), in order to represent the inspiral regime of our 
final IMR waveform model. 

Post-Newtonian perturbation theory provides a rich zoo of 
different approximants, which solve the Einstein equations up 
to a certain order of expansion, but may differ in higher order 
terms depending on how intermediate non-linear mathemat¬ 
ical operations are performed. We will use standard quasi¬ 
circular non-precessing “Taylor-approximants”, described in 
01, in particular the TaylorTl, TaylorT2 and TaylorT4 time 
domain variants, and the TaylorF2 frequency domain version. 
These are consistently derived from an energy and flux which 
incorporate up to 3.5 PN order non-spinning terms (see for 
instance 0]), 3.5 PN order in spin-orbit interaction 17^1731 . 
3PN order quadratic in spin ll74l . and 3.5PN order cubic in 
spin 1751 . 

Terms corresponding to higher PN order may also be added 
in order to conform with a particular form of the energy, flux, 
or some other property of the binary. In our work we will in 
particular use models derived from the very successful EOB 
description cma. In particular, we will use the SEOB- 
NRvl 17^ and SEOBNRv2 ||6l models in two forms, the 
proper SEOBNRvl/v2 waveform model, which is calibrated 
against numerical relativity waveforms from the SXS cata¬ 
logue ll30l . and an uncalibrated version, which we refer to as 
SEOBvl/v2, where all terms calibrated to NR have been set to 
zero. In detail, we have applied the following changes to the 
SEOBNRv2 model 161 to obtain the uncalibrated waveforms 
which we will refer to as SEOBv2: In the EOB Hamiltonian 
we have set the 4.5PN spin-orbit adjustable parametr dso 
and the 3-PN spin-spin parameter dss to zero, and we for 
non-spinning correction only keep the mass-ratio independent 
term in the parameter K, thus K - 1.712. The 3PN term 


added to the quantity 622 in the factorized waveforms is 
also set to zero, as are all NQC corrections. Furthermore, 
in order to keep the implementation of the SEOBNRv2 
in the file LALSimIMRSpinAlignedEOB. c LAE library 
0 from crashing, the parameter tStepBack was changed 
from tStepBack = 1Q8.*mTScaled to tStepBack = 
308. *mTScaled. While this work has been in progress, it 
has been found in m, that an (independent) uncalibrated 
version of SEOBNRv2 shows excellent agreement with a 
125-orbit-long mass-ratio 7 non-spinning waveform. 

D. Hybrids 

We construct hybrid waveforms, where the early time be¬ 
haviour is described by a PN approximant, and the late time 
behaviour by a NR waveform, by essentially standard proce¬ 
dures. See Refs. Il25l|60l|63]|69l for summaries of different 
methods in the literature, and Ref. Ell for a detailed descrip¬ 
tion of our specific procedure, which we briefly summarize 
here. 

First, we note that any non-precessing-binary waveform 
h{t) we consider represents an equivalence class of wave¬ 
forms, which correspond to a binary with the same masses and 
spins, and only differ by a relative time shift, a relative rota¬ 
tion, and the sense of their orbital rotation (i.e. the sign of the 
frequency w). Without restricting generality we will choose 
the frequency as positive. The PN and NR strains will 

then be related by 

/!™(f) = -H At). (2.13) 

In order to find good parameters At and (fo, we fix an initial 
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time to and fitting window length T, and minimize the quantity 

A(Af;fo,7’)= {o/^{t)-oj^^{t +At)) dt, (2.14) 

Jto 

with respect to the time shift At. Alternatively we replace ai by 
(f) to obtain a consistency and robustness check, which is par¬ 
ticularly useful when automatically processing tens of wave¬ 
forms as we do. Our original choice of oj has the advantage of 
not depending on (fQ. Other authors have replaced u) by h ll60l 
in the integrand, here we only use phase and no amplitude in¬ 
formation in order to avoid the influence of amplitude errors. 
In ITTTI we have checked that the window length T that corre¬ 
sponds to four GW cycles (two orbits) is sufficient to average 
out oscillatory waveform artefacts mostly due to residual ec¬ 
centricity. In order to hnd an appropriate phase shift (fo for 
given At from the minimization of Eq. ( |2.14| i, we align the 
phases at the middle of the matching window. 

The hybrid waveform is then defined as, 

hit) = (2.15) 


non-spinning smaller hole. Again we use SEOBNRv2 as the 
reference waveform. We see that the NR and SEOBNRv2 
waveforms agree well over the common frequency range, 
and the approximant that performs next best is the uncali¬ 
brated SEOBv2. Note that now the At time scales are much 
larger than in Eig. and that the standard PN approximants 
perform markedly worse in comparison to the NR than the 
SEOB(NR)v2 waveforms. 



- T1 

T2 

- T1-3PN 

- T4 

- SEOBvl 

- SEOBNRvI 

-SEOBv2 

- NR 


where are linear transition functions over the time win¬ 
dow t 6 [fo, to + T], 

The crucial step, however, when “glueing” PN to numerical 
waveforms is to consistently align them at some reference fre¬ 
quency, where we choose their time coordinates and phases to 
coincide. If the waveforms are identical, the relative time (and 
phase) shift between them will be independent of the matching 
frequency. We are interested in the case where the waveforms 
are different, and where we have the choice of PN approxi¬ 
mant, e.g. one discussed in Sec. |IIC| to attach to a particular 
numerical waveform. Eor non-identical waveforms the rela¬ 
tive time-shift between the PN and NR descriptions will de¬ 
pend on the frequency where the two waveforms are compared 
or matched together. A preferred PN approximant would ex¬ 
hibit a time shift relative to the NR waveform which depends 
only negligibly on time. In order to find the best-matching 
PN approximant, we compute the hybridisation time shift as 
a function of matching frequency as follows. We consider the 
time dependent frequencies upNit) and a)ffp{t). We can invert 
these functions numerically, e.g. by evaluating them for a se¬ 
ries of time steps, and then interpolating inverse functions of 
ticxj). To compute the time shift we then compute the differ¬ 
ence 


At - tp^io)) - tMRito). (2.16) 


The results are displayed in Pigs. 4]|5 Pirst, in Pig. |^we con¬ 
sider the equal mass non-spinning case, where we however 
use SEOBNRv2 and not the NR results as the reference wave¬ 
form, in order to reduce the effect of oscillations in the NR 
result due to residual eccentricity. One can see that the NR 
waveform SXS:BBH:0001 from the SXS catalogue, the un¬ 
calibrated SEOBv2 and TaylorT4 are very close to each other 
over the frequency range of the NR waveform. SEOBNRv2 
and SEOBv2 are close over the whole frequency range plot¬ 
ted, 0.01 < Ma) 22 - 


In the left panel of Pig. [^we consider the BAM mass-ratio 
18 case with a mild dimensionless spin of x\ = 0.4 paral¬ 
lel to the orbital angular momentum on the bigger BH, and a 


FIG. 4: Time shift At according to Eq. p. 16| l as a function of wave 
frequency Mo) for the case of non-spinning q = I, using SEOBNRv2 
as the reference waveform and SXS:BBH:0001 from the SXS cata¬ 
logue as the NR waveform. As is well known, the SEOBNRv2 model 
and TaylorT4 are very close to the NR result in this case. Uncali¬ 
brated SEOBv2 (uncalibrated SEOBNRv2) is also very close, while 
TaylorTl, TaylorTl cut off at 3PN order, and TaylorT2 show signifi¬ 
cant dephasing. The SEOBNRvI and SEOBvl (uncalibrated SEOB¬ 
NRvI) show rather similar performance, with good agreement with 
NR, and disagreement with SEOBNRv2 at lower frequencies. 

In the right panel of Pig. we finally consider again equal 
masses, but now with large spins of xi - X 2 = 0.98 paral¬ 
lel to the orbital angular momentum. The NR waveform is 
SXS;BBH:0172, produced with the SpEC code. The dephas¬ 
ing between the NR waveform and SEOBNRv2 is again rather 
small, although significantly larger than in the nonspinning 
case. The SEOBv2 waveform is the approximant closest to 
SEOBNRv2 for waveforms not longer than lO^M. Por lower 
frequencies, TaylorTl cut off at 3PN order is closer to SEOB- 
NRv2, but all other approximants show significantly larger 
differences. 

The time-shift analysis discussed here for three examples 
has also been carried out for a number of other cases, with 
similar results, which will be presented in more detail in a 
separate paper. The conclusion for our current work is that 
in the frequency regime relevant for hybridisation, SEOBv2 
uncalibrated EOB waveforms are almost as close to NR re¬ 
sults as the calibrated EOB model, and at low frequencies are 
typically closer to SEOBNRv2 than any other approximant at 
3.5PN order. In Ref. IfTTl we have shown a similar result, 
where we have applied a similar analysis to compare Tay¬ 
lorTl, TaylorT4 and SEOBNRvI to the BAM non-spinning 
NR waveform at mass ratio 18, and have found that SEOB¬ 
NRvI is far superior for hybridization. 

In this work we will therefore construct all our hybrid wave¬ 
forms, which form the input to our phenomenological model¬ 
ing, with the SEOBv2 model. We are thus able to calibrate an 
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FIG. 5: Left panel: Time shift Af according to Eq |2.16| as a function of wave frequency Mcj for the case of mass ratio 18 with spins {0.4,0) 
with SEOBNRv2 as the reference waveform. The PN/EOB waveforms start at Mo) = 0.01. Right panel: Same analysis for an equal mass 
binary with equal spins = 0.98. Tthe NR waveform is SXS:BBH:0172 from the SXS catalogue, see Table|^ SEOBNRvl/SEOBvl are not 
shown for this case, since the model is not valid for high spins parallel to the orbital angular momentum. 


inspiral model that is independent from the calibration process 
performed for SEOBNRv2, which will allow valuable cross 
checks. 

Frequency domain strain hybrids are constructed by sam¬ 
pling the time domain hybrid with equispaced time steps, and 
aligning all waveforms such that the value of the time coordi¬ 
nate and phase at the peak amplitude vanish. The waveforms 
are then zero-padded to all start at a time fstart and all end at 
a time fgnd, where fstart is signihcantly smaller than the begin¬ 
ning of the numerical FOB data, and fend is chosen at several 
hundred M after the amplitude peak. Different batches of hy¬ 
brids have been produced for comparisons, with the time-step 
typically chosen as 0.5M or IM. 

The time-domain amplitude for selected hybrids is shown 
in Fig.|^ The strain for the same cases is shown in Fig.|^ and 
the Fourier-domain amplitudes are shown in Figs.|^ 



t/M 

FIG. 6: Time domain amplitude of configurations at the corners of 
our parameter-space coverage, aligned at merger, showing both the 
NR waveform and the hybrid (marked as dashed lines). Thick black 
lines mark the hybridization regions. 


III. MODELING THE EINAL STATE 

While it has not been rigorously proven within general rel¬ 
ativity, the final state of the coalescence observed in all NR 
calculations is a single perturbed Kerr BH. The complex 
frequency of the “quasinotmal” ringdown provides two key 
pieces of information for modelling the GW signal, and can 
be computed directly as a function of the final mass and spin. 
As a first step to model waveforms, we will thus model the 
final state and provide simple analytical fitting formulas for 
the radiated energy and final spin as a function of symmet¬ 
ric mass ratio and an effective spin parameter. These cover a 
wider parameter space than previous work in the literature; see 
Refs. ||78l - [8Tl and references therein. This will also illustrate 
how we consfiTict analytical fits across the parameter space, 
and similar procedures will be used afterwards in our wave¬ 
form modelling. The main goal of our approach to model the 
final state is to construct accurate explicit expressions, which 
are therefore fast to evaluate and easy to generalize to larger 
calibration data sets. 

To determine the dimensionless ringdown frequency Mljro 
as a function of a/, we have interpolated the data set in 
Refs. Il82l [8^ . These well known functions are shown in 
Fig. 0 together with the wave frequency corresponding to a 
test particle on innermost stable orbit (ISCO) IIMII of Kerr 
spacetime. One can see that the steepest functional depen¬ 
dence happens for large positive spins. 

To calculate the dimensionful ringdown frequency of the 
binary, loss of energy during inspiral needs to be accounted 
for as 




(Mojrd) 

Mf 


(MoJro) 
M - ’ 


(3.1) 


and we need to model the radiated energy in addition to 
the final Kerr parameter. 

We have used final spin and radiated energy results from 
the following 106 data sets: 15 new BAM cases at mass ratios 
2, 3, 4, 8, 18, 5 old BAM cases at mass ratio 4 from pub¬ 
lished in ll85l (used only for the final spin), 50 cases from the 
SXS catalogue ||30l, and 36 cases from the RIT group Il80l . 
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FIG. 7: Real part of the time domain strain of comer cases, time coordinate is set to zero at the mid-frequency of the hybridization region, 
which is shown in red. 




Mf M f 


FIG. 8: Fourier domain amplitude of corner cases and ringdown frequency. In the left panel, the thick lines indicate the part of the NR 
waveform that was used in the final hybrid, while data indicated by a thin line contained noise or other spurious artifacts, and was discarded. 
The vertical lines indicate the peak of the amplitude, and the ringdown frequency /rq. In the right panel the leading-PN-order amplitude 
has been scaled out to illustrate detailed features of the Fourier-domain amplitude through merger and ringdown. 


For the RIT runs, two type of results are given, derived from 
the isolated horizon formalism and from the radiation, and we 
use the isolated horizon quantities due to the lower error bars 
quoted. We do not use the mass-ratio 10 results from the RIT 
data set, since the quoted result for the radiated energy dif¬ 
fers from ours and the numerical fit by a factor of two and is 
possibly a misprint. 


A. Final spin 


Before merger, the total angular momentum can be ex¬ 
pressed in terms of the individual BH spins which for our 
purposes change only negligibly during inspiral, and the or- 
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Re[W22]-Iin[(j22] . ^sco 



FIG. 9: Real (oscillatory)and imaginary (damping) part of the ring- 
down frequency as a function of the final Kerr parameter Of. A pos¬ 
itive sign denotes a final BH parallel to the initial orbital angular 
momentum, and negative sign antiparallel. 


bital angular momentum L as, 

J — L + ^ I + ^ 2 - 


(3.2) 


Due to symmetry all angular momentum vectors are orthogo¬ 
nal to the orbital plane and the only non-zero component is the 
z-component. For simplicity of notation, we therefore drop 
the vector notation, and write 7,5,, etc., for the z-component. 
We seek to approximate J, L, and correspondingly the final 
KetT parameter af - through a function of the mass ra¬ 


tio and a single effective spin parameter. Given that Eq. (3.2 1 
depends trivially on the sum of the individual spins, we at¬ 
tempt the approximation 


In the infinite mass ratio limit, rf - 0, the final spin coincides 
with the spin of the larger BH, 


7/(?7,5i, 52) = 5i = 5, af{T],x\,X2) ^Xi- 


(3.3) 


We will frequently use a rescaled version of our total spin S - 
SI + 82^0 § = 5/(1 - 2//), for which -1 < 5 < 1, consistent 
with the extreme Kerr limit. 


The data points are shown in Fig. 10 As can be seen, our 
waveform coverage is densest at equal mass and sparser at 
higher mass ratios. To produce an analytical fit, we first in¬ 
spect the data set displayed in Fig. and in particular the 
non-spinning and equal mass subsets. We find that fourth or¬ 
der polynomials in rj and S produce accurate fits for the two 
subsets, and we fix the linear term in 77 by a Taylor expansion 
around the extreme mass ratio limit as in 


df -2 y/3ri + higher order in t] and spins. (3.4) 

In order to cover the whole parameter space we extend the 
ansatz by terms quadratic in 77 and quartic in the total spin S 


to 


- foi) + S + T]'^ ^S' 
i=0,4 

f^i ^ i 


(3.5) 


!=0,4 


After fixing coefficients for consistency with the nonspinning 
and equal mass cases, the only coefficients left to fit are the 4 
numbers (fu, f 12 , fuju)- The result is, 


= S -h2V3 77 -4.399 772 4-9.397 773 - 13.181 77 ^ 
-h(-0.085 5 -h 0 . 10252 - 1 . 35552 - 0 . 8685^)77 
-h(-5.8375 - 2.09752 4 - 4.10952 4 - 2 . 0545^)772 (3.6) 


RMS errors are 6.8 x 10 2 , and 2.4 x 10 2 when comparing 
the fit with the equal spin subset. We also note that this fit 
respects the Kerr limit, i.e., < 1. 


0.25 



FIG. 10: Final Kerr parameter from Eq. |3.6| plotted as a function of 
symmetric mass ratio and total spin S. Black dots mark data points 
with equal spins, grey dots unequal spins. 


B. Radiated energy 


The data points for radiated energy are shown in Fig. 11 
together with the effective single spin fit we will now discuss. 
Inspecting the plot, clearly, a polynomial in the symmetric 
mass ratio and the effective spin will not provide the ideal 
model, and a rational function model comes to mind. Also, 
as in the final spin case, clearly a reliable accurate model of 
the whole parameter space would require further data points 
for large spins. 

It turns out that after factoring out a fit to the nonspinning 
subset, the radiated energy depends only rather weakly on the 
symmetric mass ratio. We find that the non-spinning radiated 
energy, £^^( 77 ), is very well captured by a fourth order poly¬ 
nomial, at a RMS of 3.5 x 10 2 , 

£^^( 77 ) = 0.0559745774-0.58095l772-0.960673772-H3.3524l77^ 

(3.7) 
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As can be seen in Fig. 10 the effective spin S works reason¬ 
ably well for the radiated energy. We model the dependence 
on the spin S through a simple rational function, where the 
numerator and denominator are linear in spin and quadratic 
in symmetric mass ratio, with some experimentation having 
gone into making a choice for which the nonlinear fitting pro¬ 
cedure converges well, and the result has no singularities due 
to vanishing denominator. The result of the fit, with a RMS 
is. 


error of 4 X 10 


example in Section 3 of Ref. m we obtain the SPA ampli¬ 
tude using the TaylorT4 form of the energy balance equation. 
The Fourier domain amplitude, A 22 , is then given in terms of 
the time domain amplitude A 22 and second phase derivative 
(p evaluated at the Fourier variable tf - (2nf I, where 
m - 2 for the dominant harmonic we consider. 


Mlif) - A22{tf) 



(4.1) 


l+S (-0.00303023 - 2 .OO 66 I 77 -h 1.10506jf) 

Mini ~ 1 + s (-0.67144 - 1.47569/7 + 7.30468/72) 

(3.8) 



FIG. 11: Radiated energy according to fit Eq. |3.8| plotted as a function 
of symmetric mass ratio and total spin S. Black dots mark data points 
with equal spins, grey dots unequal spins. 


IV. WAVEFORM ANATOMY AND MODEL 
A. Waveform anatomy 

7 . Amplitude 

The waveform anatomy in the time domain has been stud¬ 
ied extensively, and combined with FOB resummation tech¬ 
niques of PN results has given rise to the family of EOB- 
NR waveform models ISl [IMQl. One may distinguish 
three phases without sharp boundaries; (i) a long inspiral 
with slowly increasing amplitude, where the amplitude scales 
as M/ 7 fa/ 2/3 in tjjg low frequency limit, and the coalescence 
time as jT], (ii) an “extended-merger” characterised by a 
rapid increase in amplitude and frequency; (iii) followed by a 
damped sinusoidal ringdown. The Fourier transformation to 
obtain the frequency domain waveform can in general not be 
carried out analytically. In the low-frequency regime however, 
the stationary phase approximation (SPA) can be used to an¬ 
alytically obtain an approximate Fourier transform, which is 
used in particular to obtain the TaylorF2 PN approximant as a 
closed-form expression. Following the procedure outlined for 


In particular, to leading order the Fourier amplitude is. 


\h22\^Aof-^l\\+0{fl^)), 


Ao 



(4.2) 


In order to better emphasize the non-trivial features of the 
amplitude, we rescale our numerical data sets by the factor 
f^^jAo, to normalize all amplitudes to unity at zero frequency 
as shown in Fig. We see a structure that is sufficiently 
rich that a single analytical expression for the entire frequency 
range is difficult to achieve in terms of elementary (and thus 
computationally cheap) functions. Our strategy will thus split 
our description into an inspiral part, which models the wave¬ 
form as higher order corrections to PN expressions, a merger- 
ringdown part which builds upon the knowledge of the final 
state, and an intermediate part which describes the frequency 
regime which can not be based directly upon PN or the final 
state. 

Regarding the merger-ringdown, a crude time domain 
model, which can be Fourier-transformed analytically, is a 
sine-Exponential, which is symmetric around the peak am¬ 
plitude; 


— Tdampkl) 


fdamp ^ 


The Fourier transform \2A\ yields a Lorentzian, 


h{u) = - - 


/da: 


^{f-fRoY + h 


(4.3) 


(4.4) 


damp 


which only falls off as / 2 at large frequencies as expected 
from Eq. ( |2.8| l, due to the fact that the original time domain 
waveform hit) is only C** at the peak. Despite its oversim¬ 
plification, in particular the unphysical symmetry around the 
peak, and the incorrectly slow falloff at high frequency, the 
Lorentzian has provided a valid model for frequencies higher 
than the ringdown frequency in PhenomA/B/C models. Look¬ 
ing at Fig. 1^ the expected roughly exponential drop at high 
frequencies are clearly visible. 

A natural extension of the Lorentzian ansatz used in the pre¬ 
vious Phenom models, is to model the merger-ringdown am¬ 
plitude Amr by multiplying the Lorentzian by an exponential 
as 


^MR ^ (ys/damp) -d(f-fRD) 

Ao (/ - /rd)^ + (/dampTa)^ 


In order to find best fit parameters, we use Mathematica’s 
NonlinearModelFit function. To achieve robust conver¬ 
gence of the nonlinear least squares fit, we redefine 


4 — iy2l ifdamipy'i))-) 
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to achieve a magnitude of order unity for the new free pa¬ 
rameter 72 - The merger-ringdown amplitude model thus has 
three free parameters that are fit from numerical waveform 
data, and the real and imaginary part /RD,/damp of the QNM 
ringdown frequency are computed from the known fits for fi¬ 
nal spin and mass from Eqs. ( |3.6[ ) and p. 8 [ ) and the known 
dependence of the QNM frequencies on the dimensionless 
Kerr parameter shown in Fig|^ Interestingly, it turns out that 
73 e [1.25,1.36], i.e., it has a very narrow dynamical range 
and could also be approximated as constant 73 a; 1.3. 


The exponential factor in our ansatz Eq. (4.51 shifts the am¬ 
plitude peak from /rd for the Lorentzian to 


/max — 


/rd + 


/dampTS ^ 1 j 


72 


(4.6) 


An appropriate frequency range for fitting the model 
Eq. ( |4.5[ ) to numerical data can be determined as follows. At 
the large frequency end, we first discard frequencies with sig¬ 
nificant numerical noise by the following heuristic procedure: 
We multiply the amplitude with /^, where we find that p ^ 4 
works well. Eor physically accurate data we expect a mono¬ 
tonic decrease for the rescaled amplitude. Due to unphysical 
artefacts however, a minimum will appear, and we discard all 
frequencies higher than the one of the minimum. The result 
can be seen in Fig. While no “fine-tuning” is required, a 
good choice of frequency range for model fitting will not only 
increase the accuracy of the fit, but also lead to smoother vari¬ 
ation over the BH parameter space of mass ratio and spins, 
which will be essential when constructing a phenomenologi¬ 
cal model in Paper 2 JT]. At low frequency we find that the 
fit is only good up to the peak / = /max, while at large fre¬ 
quencies, where the amplitude is already very small, we need 
to control contamination from numerical artefacts. Also, to 
avoid boundary effects, we extend the fitting range to slightly 
smaller frequencies than / = /max- We find that two strategies 
work well: specifying the fitting range in terms of multiples 
of the ringdown frequency and the numerically determined lo¬ 
cal maximum, e.g. using the range (/max, 1-2 /rd), or using the 
width of the Lorentzian, e.g. as (/max - /damp/2, /rd -h 2/damp). 
Fits and residuals for the latter choice are shown in Figs.[l2| 
andfOl 

Since the merger-ringdown amplitude model Eq. ( |4.5| l thus 
covers a frequency range which will be covered even by the 
shortest NR simulations of BH mergers, simulations only need 
to be sufficiently long to control the systematic error of resid¬ 
ual eccentricity, to avoid low frequency artefacts from the 
Fourier transformation, and to clearly separate the unphysi¬ 
cal radiation content discussed in Sec. IIB from the physical 
part of the signal. 

Returning now to the inspiral, it is natural to model the early 
frequency behaviour as a modification of the known closed- 
form PN expression corresponding to the TaylorF2 approx- 
imant. As can be seen in Fig|^ modeling the pre-merger 
amplitude is particularly challenging for negative spins (and 
higher mass ratios), where the amplitude shows a very sharp 
drop, whereas for positive spins the frequency dependence is 



Mf 


- SXS:q1—0.95 

- SXS:q1-H-h0.980 

BAM:q8-H-H0.85 

- BAM:q18-0.8_0 

BAM:q18-^0.4_0 


FIG. 12: Comparison of ringdown fit with numerical data for 5 “cor¬ 
ner cases” as indicated in the legend. Full thin lines indicate raw 
numerical data, full thick lines indicate numerical data with a heuris¬ 
tic cutoff to remove noise at high frequencies, dashed lines indicate 
the best fit to ansatz Eq.|4.5| 



FIG. 13: Residuals between numerical data and best fits for the data 
shown in Fig.|12| 


rather shallow. In the PhenomC model 1271 the amplitude' 


taken directly in the form of Eq. (4.1 1 as resulting from the 
SPA calculation. We find however that the variation of the 
phenomenological parameters is smoother if instead of using 
the SPA amplitude defined in Eq. ( |4.1| l we re-expand the ratio 
of the separate PN expansions of A 22 (f/) and /(f/). A com¬ 
parison of different forms of the the TaylorF2 amplitude is 


displayed in Fig. 14 for the ^ = 18 case with the larger black 
hole spinning at;^^i = -0.8. One can see that the re-expanded 
forms are closer to the hybrid data. This does not hold across 
the parameter space, but it is for the negative spin cases like 
the one shown here where the difference is particularly large. 

To perform the re-expansion we set the logarithmic terms 
to zero and retain the complex terms in A 22 . Based upon com¬ 
parisons to the hybrid it was only necessary to keep terms to 
0{f^), which is equivalent to a consistent truncation at 3PN 
order. In the final result we use only the real part of this ex¬ 
pression; the contribution from the imaginary part is negligi¬ 
ble. The re-expanded SPA amplitude is given by. 


ApN(/)=Ao£l?^,■(7r/y/^ 


(4.7) 


/=0 


where Aq is the leading order / behaviour. We then cali- 
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FIG. 14: Different forms of the TaylorF2 amplitude, rescaled as in the 
right panel of Fig. are displayed for the ijr = 18 case with the larger 
black hole spinning at xi = -0.8. Shown are the numerical data 
of the hybrid waveform, the absolute value of the complex TaylorF2 
amplitude from Eq. l |4. l| l, and the real part of the TaylorF2 amplitude 
after re-expanding ton = 2.5,3,3.5 PN order. It can be seen that the 
re-expanded forms are closer to the hybrid data. Taking the real part 
or absolute value of expressions makes no significant difference for 
the case and frequency range shown. 


brate the next three natural terms in the PN expansion, 

3 

^In.s = ^PN + ^0 ^ Pi . (4.8) 

!=1 

The choice of 3 terms resulted as a compromise between accu¬ 
racy and simplicity, and adding more terms may be necessary 
in the future to further improve the low frequency behaviour, 
while further progress computing post-Newtonian terms, or 
improvements in resumming techniques may reduce the num¬ 
ber of calibration terms needed in the future. The frequency 
range for the inspiral fit was chosen by experimentation, and 
fM = 0.018 was found as a good value across our parameter 
space of hybrids, where lowering the cutoff frequency would 
not significantly increase the accuracy of fits, but increasing 
the cutoff would increase the fit residuals, as shown in Fig.flS] 
for an example case. 

Since the post-Newtonian expansion shows poor conver¬ 
gence at higher frequencies, it is not surprising that the higher 
order terms we calibrate are alternating in sign, and the prob¬ 
lem of re-constructing the signal from parameter space fits of 
the coefficients p, is poorly conditioned. We obtain signifi¬ 
cantly better results by computing the values of the fit at a 
set of 3 equispaced frequencies, and interpolating the values 
of these collocation points across the parameter space. This 
technique is used in the companion paper IJT) to construct ac¬ 
curate inspiral amplitude fits for the PhenomD model. 


2. Phase 


As already discussed regarding the construction of hybrid 
waveforms in Sec. II D[ waveforms from a system with identi¬ 
cal intrinsic physical parameters correspond to an equivalence 



FIG. 15: Residuals for the inspiral fits correspding to Eq. ( |4.8| l with 
cutoff frequencies 0.016, 0.018, 0.02. 


class of waveforms which differ by a phase and time shift. In 
the time domain, both amplitude and phase are naturally af¬ 
fected by a time shift, in the Fourier domain the amplitude is 
invariant, and the phase picks up a constant and a term lin¬ 
ear in the frequency, as can be seen by Fourier-transforming a 
shifted waveform = h{t - 


f 


V.foM = ( h(t - to)e“^‘^e dt 

^OC 

dt'. (4.9) 


_ g((0o+WZo) 


Previous phenomenological waveform constructions have 
dealt with this ambiguity by defining a standard form of the 
phase function by subtracting a fit to a linear function + 


over some suitably chosen frequency interval. In Fig. 16 


phase functions after removing different fits to Zo + w fo are 
shown, and the ringdown frequency /ro is marked. The strong 
dependence of the phase function on the alignment and the 
fact that no particular features are visible at the ringdown 
makes it difficult to learn about qualitative features of the 
phase, in particular in the merger-ringdown regime. 

For this reason we take a different path to simplify the un¬ 
derstanding of the anatomy of the phase functions, and we 
consider the first and second derivatives of the phase func¬ 
tion with respect to the frequency. For NR data, the sec¬ 
ond derivative is often noisy and we thus mainly work with 
the first derivative, which eliminates the phase shift ambigu¬ 
ity and converts the time shift ambiguity into a simple con¬ 
stant offset. In Fig. 17 we plot the derivative of the phase for 
the same 5 “corner cases” of maximal mass ratio and/or spin 
shown in Fig. and mark the ringdown frequency. We no¬ 
tice a characteristic structure rising above the background and 
centered at the ringdown frequency, which reminds us of the 
Lorentzian Eq. ( |4.4| i we have previously employed to model 
the amplitude function. Taking a derivative of the Lorentzian, 
we can derive that the ma xima of its derivative are located at 
fRD ± /damp I V3- In Fig. 


IVA2 


we compare this prediction 
for the second phase derivative with NR data from a BAM 
simulation, where we can compute a relatively clean second 
phase derivative, and we find excellent agreement. In order 
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Mf 

FIG. 16: Fourier domain phase corresponding to a hybrid constructed 
for SXS NR data corresponding to equal mass data with equal spins 
of Xi = 0.98, with a term (pQ + a) Iq subtracted over a frequency 
range. The frequency ranges of the fits are indicated in the plot, 
where Mfs = 0.0035 is the start frequency of clean hybrid hybrid 
data, and M/mr - 0.018 the frequency where we typically separate 
between inspiral and merger-ringdown. 

to capture the “background” in which the Lorentzian is em¬ 
bedded, we add terms with negative powers of the frequency 
/, to capture the steep inspiral gradient, and our basic model 
becomes. 


frequency range with a single or at least only 2 fitting regions. 

For consistency with our amplitude model, we again choose 
Mf - 0.018 as the transition between our inspiral and high- 
frequency descriptions. For the merger-ringdown part, we find 
that a simple choice in terms of the ringdown frequency, e.g. 
[0.45,1.15] /rd is sufficient. The upper frequency 1.15/rd ap¬ 
proximates the highest frequency for which we have clean NR 
data. For the full IMR phase we use the above fit for frequen¬ 
cies larger than 0.5 /rd. 

For the intermediate phase we then choose a fitting window 
of [0.018,0.6/rd], and we use 2 fitting coefficients p 2 - 1, 
P3 = 4. 



0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 

Mf 


/mr " O'! + X ^ 

i=2 


fLn+if-fRor 


(4.10) 


where the constant ao acts as an overall time shift of the phase 
function, which could be determined by matching to a de¬ 
scription of the inspiral, and the a, and a are free parameters 
which can be fit to numerical data. The number of free param¬ 
eters a,- and the choice of powers p„ defines a specific merger- 
ringdown ansatz. In order to find a good choice for the number 
of parameters and we first produce fits with the free param¬ 
eters (a'i,Q' 2 ,p 2 ), and determine the best fit power /?„ across 
the parameter space. We find that p 2 ~ 2, so we choose the in¬ 
teger p 2 -2 as, the leading contribution. Further terms can be 
added to the ansatz, e.g. until the residual across the parameter 
space resembles pure noise, and no further information can be 
extracted via fits. The specific choice made for the PhenomD 
model discussed in the companion paper is to add one more 
term with p 3 = -1/4. In order to improve the convergence 
of Mathematica’s NonlinearModelFit function, we replace 


/damp by aofdamp in Eq. (4.10l and also fit to ao- We find that 
the parameter ao is in the range [0.98,1.04] and compensates 
for errors in our final mass and final spin fits near the bound¬ 
ary of our calibration region. For the actual PhenomD model 
discussed in IT] this parameter will be set to unity. 

For intermediate frequencies, roughly between Mf - 0.018 
and half the ringdown frequency one finds that the best choice 
of p 2 for a single power law coefficient in Eq.(4.10l \s p 2 ~ 1, 


so we choose the integer p 2 - I as the leading contribution 
in this regime, and we break up our model into different fre¬ 
quency ranges in each of which we use a very simple ansatz. 
A more sophisticated ansatz may successfully treat the entire 


FIG. 17: Phase derivative for “corner cases” as indicated in the 
legend, vertical straight lines mark the ringdown frequency fpo. 



FIG. 18: Second phase derivative for the BAM waveform with q = 2, 
Xi = 0.75, vertical straight lines mark the ringdown frequency fpo 
and fpD ± fdamp I V 3 - 


We thus have succeeded in directly using information about 
the final spin and mass in the ringdown description of both 
amplitude and phase, and turn now to using PN information 
for the inspiral. We can treat the inspiral along the lines of the 
amplitude in Sec jlV AT For the amplitude we have added the 
next three higher order PN terms to the the known TaylorF2 
expression. For the phase, the 4PN terms corresponds to a 
constant offset for the phase, to which we add the next three 
non-trivial higher order PN terms, thus four PN terms to be 
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fitted in total. 


V. SUMMARY 

In this paper we have discussed preliminary work that will 
enable us to construct a new phenomenological waveform 
model for non-precessing waveforms in a companion paper 

m. 

We have first presented new NR simulations which include 
data sets with a range of spin values at mass ratios 4, 8 and 18, 
which allow us to signihcantly extend the calibration range 
of current non-precessing inspiral-merger-ringdown models. 
We have demonstrated sixth order convergence for one of the 
^ = 18 cases and computed a phase error estimate on the order 
of 0.05 radians after aligning near merger. We hnd in particu¬ 
lar small secular drift, and errors that are dominated by noise 
and oscillations. We have illustrated high accuracy also for a 
~10 orbit aligned-spin case, where we have not completed the 
convergence series to save computational cost, and conclude 
that finite difference codes based on the “moving-puncture” 
paradigm ll38l provide a very robust tool to obtain NR 
data sets in new regions of parameter space. 

We have then discussed our procedure to construct hybrid 
PN/numerical relativity waveforms in the time domain, and 
how to convert them to the frequency domain. Special em¬ 
phasis has been put on the alignment procedure which is at the 
core of the hybridization procedure, and we have presented a 
comparison of our NR data and different PN approximants, 
from which we conclude that the most appropriate approx- 
imant to use in our hybrid construction is the uncalibrated 
SEOBv2 model. We have also described our simple criteria to 
determine an appropriate frequency range for the htting pro¬ 
cedure expressed in Eq. ( |2.10| i, which determines the optimal 
alignment of PN and NR waveforms. The PN comparison 
has been based on a time-shift analysis directly relevant to the 
hybridization procedure, which we have however also found 
very useful for other applications, such as studying the time- 
dependent dephasing of numerical simulations due to numer¬ 
ical error. 

A crucial part of a waveform model is a prediction of the 
final mass and spin from the initial configuration, which then 
allows one to compute the complex quasi-normal ringdown 
frequencies. While there are a number of numerical hts to h- 
nal mass and spin available in the literature ll78l - l8Tl . we have 
presented new hts here, which extend the calibrated range, 
are particularly simple, and can easily be extended to new 
data sets or higher accuracy as further improvements of our 
waveform model are developed. An essential feature of the 
hnal spin ht for waveform modeling is that it respects the ex¬ 
treme Kerr limit. Also, the hnal state hts illustrate techniques 
for constructing parameter space hts which we will use in 
the companion paper HI to construct the PhenomD waveform 


model. 

Previous phenomenological waveform models have used a 
Lorentzian to model the ringdown part of the Eourier am¬ 
plitude, which however leads to an incorrect high frequency 
falloff behaviour instead of (9(/“^)VM). Here we show 
that not only does a simple exponential factor lead to an ex¬ 
cellent agreement with numerical data, but we also identify a 
Lorentzian of the same width, centered at the ringdown fre¬ 
quency in the derivative of the frequency domain phase. 

Eor low frequencies, we model the amplitude and phase by 
the standard TaylorE2 approximant, augmented by pseudo-PN 
terms tuned to our hybrids constructed from SEOBv2 and NR 
data. Eor the amplitude, we use a re-expanded form with the 
imaginary part of the amplitude set to zero, i.e. we neglect the 
very small deviations between the GW phase and twice the 
orbital phase some of us discussed in Ref. ifTTl . We determine 
an appropriate cutoff frequency, where our inspiral model is 
separated from our merger ringdown model as Mf - 0.018, 
considering the dependency of the inspiral amplitude ht resid¬ 
ual on this cutoff frequency. 

In order to model how the phase Lorentzian is embedded in 
the background functional dependance of the phase, we pro¬ 
pose to complete the model by adding inverse powers of the 
frequency and identify the most appropriate leading term. 

Having identihed the Lorentzian feature in the phase, and a 
frequency up to which even a simple extension of the TaylorE2 
approximant yields a good ht to the data, paves the way for 
complete waveform models of extremely high hdelity to the 
physical signal. 

A hrst incarnation of our approach to develop phenomeno¬ 
logical waveform models that are sufficiently accurate for 
the advanced detector era is the PhenomD model, which we 
present in a companion paper ifTl . 
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